Study on Bingham fractional damage model of backfill material under different moisture content conditions

Filling mining technology is an important representative technology to realize green and low-carbon mining. The backfill materials have distinct rheological characteristics under the long-term action of formation loads and groundwater seepage. In order to study the creep characteristics of backfill materials under different moisture contents and reveal their aging-mechanical properties, based on the Riemann-Liouville fractional calculus and damage mechanics theory, the fractional element and damage variables are introduced to improve the traditional Bingham model, and the fractional Bingham creep damage model is proposed. Based on the experimental data of gangue cemented backfill under different moisture content, the parameters of the creep model are obtained by using user-defined function fitting and the least square method. The results show that the improved Bingham fractional creep damage model can describe the whole creep process of backfill materials under different moisture contents, and the rationality of the model is verified. Compared with the traditional Bingham model, the fitting degree of the Bingham fractional creep damage model is higher, which solves the problem that the traditional Bingham model cannot describe the nonlinear creep stage. Model parameter α and ξ increase with the increase of axial stress and moisture content. Under the same moisture content, η gradually increases with the increase of axial stress. This work has a certain reference significance for studying the mechanical properties and creep constitutive model of backfill materials containing water.


Introduction
Coal is an important basic energy for Chinese industrial economy and social development.It can be seen from Fig 1 that annual global coal consumption quantity is above 150EJ from 2011 to 2022.And global coal consumption increased by 9.95% in the year 2020 to 2021 [1].
On March 17, 2022, the National Energy Board issued the "2022 energy work guidance" which pointed out that it is necessary to consolidate the foundation of energy supply security, accelerate the green and low-carbon transformation of energy, and enhance the flexibility and resilience of the energy supply chain.Therefore, coal green and low-carbon development is imperative.Coal-based solid waste-filling mining technology is an important representative technology to realize green and low-carbon mining.It has remarkable technical advantages in surface subsidence control, ecological environment protection, mine solid waste disposal and utilization, and green and low-carbon emission reduction [2].
Filling mining technology has been widely used in mining engineering.The roof is supported by backfill materials, which greatly reduces the roof deformation and the ground settlement, and has a good control effect [3].However, the backfill materials has obvious creep characteristics under the long-term effect of stratum load and groundwater seepage.Therefore, it is of great engineering significance to study the creep characteristics of backfill materials under different moisture contents and reveal its time effect-mechanical properties for the evaluation of long-term stability of backfill materials.
At present, many scholars have done a lot of researches on the mechanical properties of the different backfill materials, and have achieved plenteous research achievements.Through the seepage test of coal gangue materials, Li [4] found that that adding a certain concentration of coal fly ash can achieve a better water barrier.Nie [5] revealed that significant differences in the mechanical properties, acoustic emission (AE) responses, maximum principal strain fields, and ultimate failure modes of CGBS under different static-dynamic combined loading paths.Brett [6] establishes the link among induced curing pressure, multiphysics processes, and interface behavior for FRCC materials through analysis of the state-of-the-art research findings on the FM-ITZ of FRCC materials.Gao [7] utilized industrial graphene oxide (GO) and fly ash (FA) to create a high-performance, low-cost, and environment-friendly material for backfill.Sadat [8] discussed the possibility of utilizing the construction and demolition waste (CDW) as cemented recycled aggregate backfill (CRAB) materials in the underground mining industry.Yang [9] revealed the influence of curing age on the dynamic mechanical properties and deformation and failure characteristics of backfill materials under intermediate strain rate.Ran [10] used acoustic emission technology to monitor the damage evolution of strip and column cemented gangue backfill bodies.Rafat [11,12] found that the tremendous anti-blast performance of the tubular column filled with concrete commonly known as the CFST column as a wise substitute for the conventional RC column.Zhu [13] studied the influence of different mixing ratios on the rheological properties of cemented gangue backfill, and determined the appropriate dosage range.However, the above research did not consider the effect of moisture content on the creep characteristics of backfill materials.At present, the research on the effect of moisture content on the creep properties of materials is mostly focused on rock materials.Li [14] deeply analyzed the relationship between coal (rock) permeability and the moisture content through seepage experiments of coal (rock) with different moisture contents.Duan [15] used a heatfluid-solid triaxial servo seepage device to measure the permeability of coal (rock) under different moisture contents and revealed the effects of the slippage effect, permeability, and moisture content.Chen [16] analyzed the acoustic emission damage characteristics of coal (rock) under hydraulic action through uniaxial compression tests under different moisture contents.Wang, [17,18] established a creep damage model of water-bearing oil shale.Triaxial creep tests of oil shale at different moisture contents were conducted.The rationality of the model was then verified.According to the variation relationship of sandstone strength, the relative strength criterion of sandstone under hydro-mechanical coupling was established [19].
The backfill materials material is different from the rock material, and their creep characteristics are obviously different, and the creep of backfill materials often shows a strong time effect.At present, there have been few studies on the creep constitutive models of backfill materials containing water.It is of great significance to establish a creep model that can describe the whole creep process of backfill materials containing water.Various creep models, stress; x0Dx, the Riemann-Liouville fractionalorder integral operator; εije, the strain of the elastic element; εijf, the strain of the fractional element; εijvp, the strain of the viscous element with a damage variable; εijv, the strain of the viscous element; εijt, the total strain; a and b, the material parameters; c, the cohesive force of the backfill materials; E, the elastic modulus; e ij , the deviatoric strain tensor; F, the backfill materials yield function; G 1 , the shear modulus; I 1 , the first invariant of the stress tensor; J 2 , the stress deviator second invariant; K 1 , the bulk modulus; n, the inherent material coefficient; Q, the plastic potential function; S ij , the deviatoric stress tensor; t m , the creep failure time of the object; α, the order number of the fractional calculus; Γ, the Gamma function; δ ij , the Kronecker tensor; ε, the One-dimensional axial strain; ε ij , the strain tensor; ε m , the average strain in three directions; η, the viscosity coefficient; θ, the internal friction angle of the backfill materials; μ, the damaged element viscosity coefficient; ξ, the the inherent coefficient of the fractional element; σ, the One-dimensional axial stress; σ ij , the stress tensor; σ m , the average stress in three directions; σ s , the yield strength of an object; ψ, the constant.
including the Maxwell model [20], Kelvin model [21], Bingham model [22], Burgers model [23] and Nishihara model [24], are currently widely used.The Bingham model includes viscous, elastic, and plastic elements, The Bingham model can well describe the steady creep stage of the backfill materials.Compared with these four classical models, the Bingham model has the advantage of fewer parameters and can describe the instantaneous deformation and steady creep well.However, owing to the limitations of its constitutive equation, it can only describe the linear creep stage and cannot accurately describe the nonlinear creep stage.
A large number of creep phenomena show that the strain state of a point in the material is not only related to the same instantaneous stress state of the point, but also to the whole stress history of the point before this time.By coincidence, the fractional time derivative is actually a differential-integral convolution operator, and the integral term in its definition fully reflects the historical dependence of the development of the system function and is a powerful mathematical tool for modeling processes with strong memory [25].In recent years, a fractional element based on Fractional Order Theory has been gradually applied to the study of creep models.Zhang [26] proposed a triaxial creep model of deep coal based on fractional derivative considering the temperature effect for triaxial stress state conditions.Based on the modified Nishihara model, Cheng [27] proposed a nonlinear creep model, considering the viscoelastic-plastic and damage evolution characteristics of the rock.On the basis of the Caputo variable-order fractional derivative, Liu [28] proposed the Caputo variable-order fractional creep model.In order to determine the nonlinear creep characteristics of rock under cyclic loading and unloading conditions, Zhang [29] proposed the nonlinear Kelvin model and damage viscoplastic model.Based on the fractional order theory, Huang [30] established a viscoelastic-plastic model considering the initial damage of coal samples.Gao [31] established a new variable fractional rheological model to describe the full-stage creep behavior of rock.
Based on the above problems, the fractional element [32] proposed previously was introduced in this paper.Its constitutive model has fewer parameters than other models and allows for convenient calculations.The strain of the fractional element has a power function relationship with time, which is consistent with the initial creep stage in the creep process of backfill materials.The model reflects the nonlinear characteristics of the creep phenomenon.To characterize the accelerated creep stage, a viscous element with damage variables was introduced.The above elements were combined with the traditional Bingham model to form an improved Bingham fractional creep model.One-and three-dimensional constitutive equations of the model under three different conditions ( _ εðt !1Þ ¼ 0; s < s s , _ εðt !1Þ > 0; s < s s , _ εðt !1Þ > 0; s � s s ) were deduced, and a simple parameter identification method was created.Based on the experimental data of gangue cemented backfill under different moisture content [33], the parameters of the triaxial creep data under different moisture contents were then identified, and the test curve was compared with the theoretical curve to verify the rationality and effectiveness of the model.Correlations between the parameters of the model and the moisture content were analyzed.The research results can provide an important reference for the study of the influence of moisture content on the mechanical properties of backfill materials and the creep constitutive model.

Fractional calculus
Fractional calculus has evolved on the basis of integer order calculus, evolving many definitions such as Riemann-Liouville, Captuo and Grunwald-Letnikov.Riemann-Liouville fractional calculus [34] is the earliest defined and most complete fractional calculus.The most widely used of these is the Riemann-Liouville fractional calculus.Moreover, the definition of Riemann-Liouville fractional calculus is more rigorous in mathematical expression.To facilitate the model derivation and numerical calculations, Riemann-Liouville fractional calculus was used in this study.The integral of order α of a function f(x) is defined as follows: The fractional derivative is defined as follows: Where the lower left and lower right indices of x 0 D x indicate the ranges of integration, α is the order number of the fractional calculus (0 < α, m − 1 < α < m, m 2 N*), and Γ is the Gamma function, where GðzÞ ¼ The Laplace transform formulas for fractional calculus are as follows: Where the Laplace transform of f(t) is denoted as � f ðiÞ.

Fractional element
The state of an object is assumed to be between an ideal solid and an ideal fluid.The relationship between stress and strain is expressed by a fractional derivative, and this object is referred to as a fractional element [32].The mechanical model is shown in Fig 2.
The constitutive equation of the fractional element is as follows [35,36]: Where ξ is the inherent coefficient of the fractional element.When α = 1 and sðtÞ ¼ x dεðtÞ dt , dεðtÞ dt is the rate of strain, and the fractional element is equivalent to a damper element, which represents an ideal fluid.When α = 0 and σ(t) = ξε(t), the fractional element is equivalent to a spring element, which represents an ideal solid.
When σ(t) = const, the fractional element describes the creep behavior under the condition of constant stress.Eq (4) is integrated according to Riemann-Liouville fractional calculus theory, and the creep equation of the fractional element can be calculated as follows: From this, it can be determined that the strain of the fractional element has a power function relationship with time.The creep curves for different α values are plotted according to Eq (5) in Fig 3, showing that the fractional element can characterize the initial creep stage.
Similarly, when ε(t) = const, a fractional element will describe the stress relaxation of the creep motion.At this time, the stress of the fractional element changes with time, and the strain remains unchanged.The relaxation equation of the fractional element can be derived as follows:

Viscous element with damage variables
When the backfill materials containing water enter the stage of accelerated creep, internal damage will appear and gradually develop into macroscopic cracks.To accurately describe the The damage variable first appeared in damage mechanics.The evolution equation of the damage variable proposed by the former Soviet Union scholar Kachanvo is [37]: The creep failure time can be obtained by integrating the Eq (7): Combining Eqs ( 7) and ( 8), the evolution equation of damage variable with time is: Therefore, the constitutive equation of the viscous element with damage variables is as follows: Where σ s is the yield strength of an object, σ is the applied stress, μ is the damaged element viscosity coefficient, n is the inherent material coefficient, and t m is the creep failure time of the object.
The creep curves for different values of n when s ¼ 60MPa; s s ¼ 50MPa;

Bingham model
The traditional Bingham model [22] is composed of a Hooke element and an elastic viscoplastic element in series.The mechanical model is shown in Fig 6.
The constitutive equations are as follows: Where ε is the strain, _ s and _ ε are the rates of stress and strain, respectively, E is the elastic modulus, and η is the viscosity coefficient.the applied stress does not reach the yield strength of the backfill materials, the stable creep stage is entered.The strain increases linearly with time, and the backfill materials are viscous.This stage is described by the viscous element ③.If the applied stress reaches the yield strength of the backfill materials, the accelerated creep stage is entered.There is a nonlinear relationship between the strain and time.The overall trend is concave downward, and damage occurs inside the backfill materials mass.This stage is described by a viscous element containing damage variables.

One-dimensional creep equation of improved Bingham model
The improved creep model satisfies the following conditions: (1) When _ ε t ! 1 ð Þ ¼ 0; s < s s , the model is equivalent to a Hooke element and a fractional element connected in series.The stress of the backfill materials has not reached its yield strength, and the deformation is in the initial creep stage.According to the principle of elements in series, the stress of each element is the same, and the total strain of the model is equal to the sum of the strain of each element.The following constitutive equation of the improved creep model can be obtained: Where σ is the initial stress, and E is the elastic modulus.
(2) When _ ε t ! 1 ð Þ > 0; s < s s , the model is equivalent to a Hooke element, fractional element, and viscous element connected in series.The creep reaches the stage of steady creep.The constitutive equation of the creep model is all parts of the model are involved in the creep.According to the superposition principle, the constitutive equation of the improved creep model is as follows: By combining Eqs ( 13)-( 15), the fractional viscoelastic-plastic one-dimensional creep model expression is obtained as follows:

Three-dimensional creep equation of improved Bingham model
In actual working conditions, such as backfill materials mining and excavation, backfill materials are mainly in a three-dimensional stress state.To allow the model to have a certain practical reference significance in practical engineering, the creep equation under the three-dimensional stress state is now deduced.The theoretical basis for previous work is deduced [38], and the following assumptions are made: 1. backfill materials is an isotropic material with consistent damage in all directions.
2. Damage occurs only in the accelerated creep stage and does not occur in other stages of creep.
3. The damage duration is consistent with the corresponding creep time.
Under the three-dimensional stress state, it is assumed that the total strain produced by the creep of backfill materials is ε t ij .The strain of the elastic element is ε e ij , the strain of the fractional element is ε f ij , the strain of the viscous element is ε v ij , and the strain of the viscous element with a damage variable is ε vp ij .According to the superposition principle, According to elastic mechanics, the elastic element has the following relationship.The stress tensor σ ij can be divided into a spherical stress tensor σ m δ ij and a deviatoric stress tensor S ij .σ m is the average stress in three directions.The strain tensor ε ij can be divided into a spherical strain tensor ε m δ ij and a deviatoric strain tensor e ij .ε m is the average strain in three directions.Thus, the following relationship holds: where δ ij is the Kronecker tensor [39], which is According to the generalized Hooke's law, where K 1 is the bulk modulus, and G 1 is the shear modulus.
Therefore, the creep equation of an elastic element under a three-dimensional stress state is The spherical stress tensor has little effect on creep, ignoring the effect of spherical stress tensor on creep.The creep equation of a fractional element under a three-dimensional stress state is The creep equation of a viscous element under a three-dimensional stress state is The creep equation of the viscous element with damage variables under the three-dimensional stress state is Where , F is the backfill materials yield function, F 0 is the initial value of the backfill materials yield function, which is usually set to 1, a constant, which is usually set to 1, and Q is the plastic potential function.According to the associated flow law, Q = F.According to the superposition principle, the three-dimensional creep equation is obtained as follows: Based on the viscoelastic-plastic model, the Drucker-Prager yield criterion [40] is adopted, and its equation is Where J 2 is the stress deviator second invariant, I 1 is the first invariant of the stress tensor, and a and b are material parameters.a and b are defined as follows: a ¼ sin y ffi ffi ffi 3 p ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Where θ is the internal friction angle of the backfill materials, and c is the cohesive force of the backfill materials.
In the three-dimensional axially creep test of backfill materials, the following relationships are satisfied: By combining Eqs ( 17)-( 27), according to the superposition principle, the viscoelastic-plastic damage axial creep equation of the backfill materials containing water under the threedimensional stress state is obtained as follows: Parameter identification and model validation

Parameter identification
(1) Determination of G 1 and K 1 In Eq (28), The elastic modulus E and Poisson's ratio ν were measured by uniaxial compression tests, and G 1 and K 1 were obtained by substituting into the above equation.
(2) Determination of ξ and α The data in the initial creep stage were fitted with the following nonlinear functions: The parameters of the model were obtained simultaneously as follows: The Γ function value was calculated by the Maple software [41].
(3) Determination of η The data in the steady creep stage were fitted with linear functions: The parameters of the model were obtained simultaneously, and the following formula was applied: (4) Determination of t m , μ, and n t m was obtained from the creep test, and the parameters μ and n are obtained by the following method.It was assumed that the functional relationship between the strain and parameters was ε i = f(t i , μ, n).There were a total of m groups of experimental data (t k , ε k )(k = 1,2,3,� � �,m).First, a set of initial values of parameters μ 0 and n 0 was specified, so that the initial parameters of the model were determined.Each value of t k was substituted into the creep equation to obtain the theoretical value εk .To minimize the error, the calculation method reported previously was adopted [42].When 2 was the smallest, the following conditions must be met: When using this method, iterative divergence will occur if the initial value is not selected properly.Therefore, a "damping factor" d was added to ensure the method iteratively converged [38].

Model validation
To deeply explore the effect of water on the mechanical properties of backfill materials, triaxial creep experiments of backfill materials under different moisture contents were conducted.The experimental data came from the triaxial creep test of gangue cemented filling material under different moisture contents.The detailed experimental procedure was reported previously [33].
According to the test procedure of the International Society for Rock Mechanics (ISRM), a standard cylinder sample with a diameter of 50 mm and a height-to-diameter ratio of 2:1 was prepared.To analyze the influence of the moisture content changes on the mechanical properties of the backfill materials, four groups of samples with moisture contents of 0% (dry), 8.8%, 17.6%, and 22% (completely saturated) were prepared by using a vacuum saturator and a drying oven, respectively.There were 24 samples in each group of six.The long-term strengths of the samples with 0%, 8.8%, 17.6%, and 22% moisture contents were 7.19, 6.22, 4.62 and 4.14 MPa, respectively.The composition of gangue cemented filling material is shown in Table 1.The basic physical and mechanical parameters of the samples under different moisture contents are shown in Table 2.
The MTS 815.02 electro-hydraulic servo-controlled rock mechanics testing system was used in the test equipment, and the maximum axial compression can reach 1700 kN.The compressive axial load rate can be controlled within 10 −5 ~1 mm/s.The equipment can carry out creep tests.
In this test, the confining pressure was first loaded to 1 MPa, and then the confining pressure was fixed.Axially controlled loading was then applied.The loading rate was 0.02 kN/s, the incremental load was 4 kN (2.04 MPa) for each stage, and the next stage was loaded every   As can be seen in Fig 9, with the increase of axial stress, the total deformation of the specimen increased gradually under the condition of a certain moisture content.Under the action of low axial stress, the specimen only had immediate deformation and initial creep stage.With the increase of axial stress, the specimen gradually underwent a constant creep stage and accelerated creep stage.It showed that the creep rate increased with the increase of axial stress.Under certain axial stress, with the increase of moisture content, the total deformation of the specimen increased gradually, and the total creep time decreased gradually.It showed that the moisture content had an obvious deterioration effect on the specimen, which greatly accelerated the creep failure process of the specimen.
According to the experimental results, using the parameter inversion method, the inversion results of the creep parameters are shown in Table 3.When the confining pressure was 1 MPa and the moisture content was 0%, the constant-velocity creep rate of backfill materials was basically 0 at the first six stress levels.At this time, the first formula of Eq (28) was used.Under the 7th and 8th stress levels, the backfill materials strain included not only elastic and viscoelastic strains but also viscous strains.The second formula of Eq (28) was used.The long-term strength of backfill materials exceeded the 9th stress level.The third formula of Eq ( 28

Relationship between parameters of model and moisture content
According to the data of Table 2, the relationship between parameters ξ, α, η and moisture content and axial stress is drawn, as shown in

Fig 5 .
Fig 5.As shown in Fig 5, the viscous element with damage variables can theoretically characterize the accelerated creep stage.

Fig 4 .
Fig 4. Viscoplastic element with damage variables.https://doi.org/10.1371/journal.pone.0295254.g004 To overcome the shortcomings of the traditional Bingham model, the fractional element is introduced.The fractional element and viscous element are connected in series with the Bingham model.Then, the viscous element in the Bingham model is replaced by a viscous element with damage variables.An improved five-element Bingham model is proposed.As shown in Fig 7, element ①, element ②, element ③, element ④, and element ⑤ are arranged from left to right.The creep curve is shown in Fig 8.The whole process of creep is divided into three stages: initial creep, steady creep, and accelerated creep.Different elements play their respective roles in different creep stages.The instantaneous strain of the backfill materials mass is described by elastic element ①.Entering the initial creep stage, the relationship between strain and time is nonlinear, and the overall trend is upward convex.The fractional element strain and time followed a power function relationship, and this stage is described by the fractional element ②.If

Fig 8 .
Fig 8. Classical creep curve.https://doi.org/10.1371/journal.pone.0295254.g008 7200 s until the specimen failed.After the specimen was destroyed, the confining pressure was first unloaded to zero, then the axial load was unloaded to zero, and finally, the test data was derived.Fig9shows a graph of the results of the triaxial creep test performed on the sample.Accelerated creep occurred when the moisture content of the sample was 0% and the axial stress was 18.36 MPa, and the moisture content was 8.8% and the axial stress was 14.28 MPa.The moisture content of 17.8% and 22% did not show the accelerated creep stage, because the interval of axial stress was large, and the failure occurred before reaching the next axial stress.

Fig 10 .
Fig 10.Fitted creep curves of backfill materials with different moisture contents [33].(a)0.00%(b)8.8%(c)17.6%(d)22.0%.https://doi.org/10.1371/journal.pone.0295254.g010 Fig 11.It can be seen from Fig11that under the same axial stress, α increases gradually with the increase of moisture content.Under the same moisture content, α gradually increases with the increase of axial stress.In general, α increases with the increase of axial stress and moisture content.Under the same axial stress, ξ increases gradually with the increase of moisture content.Under the same moisture content, ξ gradually increases with the increase of axial stress.In general, ξ increases with the increase of axial stress and moisture content.Under the same moisture content, η gradually increases with the increase of axial stress.